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Abstract 

Traditional chemical kinetics may be inappropriate to describe chemical reactions in 
micro-domains involving only a small number of substrate and reactant molecules. 
- - ■ Starting with the stochastic dynamics of the molecules, we derive a master-diffusion 

f-*) i equation for the joint probability density of a mobile reactant and the number of bound 

' substrate in a confined domain. We use the equation to calculate the fluctuations in the 

' number of bound substrate molecules as a function of initial reactant distribution. A 

^ , second model is presented based on a Markov description of the binding and unbinding 

Q I and on the mean first passage time of a molecule to a small portion of the boundary. 

These models can be used for the description of noise due to gating of ionic channels 
. by random binding and unbinding of ligands in biological sensor cells, such as olfactory 

I cilia, photo-receptors, hair cells in the cochlea. 

> ■ 

g ; 1 Introduction 

(N I 

T^ij- ' Biological micro-structures such as synapses, dendritic spines, subcellular domains, sensor 

O ■ cells and many other structures, are regulated by chemical reactions that involve only a 

^ . small number of molecules, that is, between a few and up to thousands of molecules. A 

^! chemical reaction that involves only 10 to 100 proteins can cause a qualitative transition in 

^ [ the physiological behavior of a given part of a cell. Large fluctuations should be expected 

^ ' in a reaction if so few molecules are involved, both in transient and persistent binding and 

H ■ unbinding reactions. In the latter case large fluctuations in the number of bound molecules 

'i> • should force the physiological state to change all the time, unless there is a specific mechanism 

^ , that prevents the switch and stabilizes the physiological state. Therefore, a theory of 

H ' chemical kinetics of such reactions is needed to predict the threshold at which switches 

■ ■ ' occur and to explain how the physiological function is regulated in molecular terms at a 

sub-cellular level. 

A physiological threshold can be determined in molecular terms, for example, when the 
number of activated molecules exceeds a certain value. The standard theory of chemical 
kinetics is insufficient for the determination of the threshold value, because it is based on the 
assumption that there is a sufficiently large number of reactant molecules and it describes 
the time evolution of only the average number of molecules. The standard theory of reaction- 
diffusion describes chemical reactions in terms of concentrations so that fluctuations due to 
a small number of molecules are lost. 
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For example, in dendritic spines of neurons a flow of calcium entering through the NMDA 
channels can induce a cascade of chemical reactions. As calcium ions diffuse they can bind, 
unbind, and leave the spine without binding. But if enough calcium binds specific molecules, 
such as calmodulin, then certain proteins become activated, such as CaMK-II, which are 
involved in regulating synaptic plasticity . Now if sufficiently many of them are activated 
at about the same time and thus the threshold is exceeded, additional changes can be induced 
at the synapse level, affecting the physiological properties of a neuron. In particular, such 
changes may include a modification of the biophysical properties of some receptors and/or 
increase the number of channels at a specific area of the synapse, called the postsynaptic 
density. It is unclear how many CAMK-II are needed for crossing the threshold, but the 
range is somewhere between 5 to 50. It is remarkable that as few as 5-50 molecules can 
control the synaptic weight The number of phosphorylated CAMK-II, activated after a 
transient calcium flow, depends on the location of the proteins, the location of the channels, 
the geometrical restrictions imposed by the spine shape, and the state of the proteins. All 
of these factors play a crucial role in determining the threshold. 

The photoreceptor cells are another example, where fluctuations in the number of bound 
molecules determine the physiological limitations of the cell. Indeed, in the outer segment of 
cones and rods of the retina, the total number of open channels fluctuates continuously due to 
the binding and unbinding of specific gating molecules to their receptors. These fluctuations 
are directly converted into fluctuations of the membrane potential, which are called "dark 
noise" , and thus determine the signal to noise ratio |2j for a photon detection. The fluctuation 
in the number of open channels is regulated by the number of gated molecules and depends on 
the geometry of the outer segment and the distribution of channels. It is not clear what are 
the details of the biochemical processes involved in regulating the number of open channels, 
but interestingly, the signal due to a single photon event is sufficient to overcome the noise 
amplitude in rods, but not in cones, although their biochemical properties are similar. As 
the binding molecules diffuse in the cell they can bind and unbind to channels, thus causing 
fluctuations in the number of open channels. The fluctuation depends on the arrival time of 
the binding molecules, called cGMP, to the channel binding sites. 

In the mathematical description of the binding and unbinding reactions, we model the 
diffusion of the particles as Brownian motion and binding occurs when a particle reaches a 
binding site. The binding probability depends on the geometry of the domain and on the 
distribution of the channels. A channel opens when it binds to two or three cGMP molecules 
and in the absence of light, the number of open channels is small, approximatively 6 to 10 
per micro-domain in a mammalian rod, when there are only 60 cGMP molecules. Due to 
the random binding and unbinding of the molecules to the channels the number of open 
channels is a stochastic process. 

In this paper we start with the stochastic dynamics of the reactant molecules in a micro- 
domain and derive a master-diffusion equation for the joint probability density of the mobile 
reactant and the number of bound immobile substrate molecules. We use the equation 
to calculate the fluctuations in the number of bound substrate molecules as a function of 
initial reactant concentration. We apply the present theory to the computation of the mean 
and variance of the fluctuation in the number of open channels, and find their dependence 
on the initial number of the mobile reactant, the geometry of the compartment, and the 
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distribution of channels. 

Our model can predict the fluctuation intensity produced by binding and unbinding to 
channels in conflned compartments, such as the a compartment (the space between two con- 
secutive disks) of cone and rod outer segment, or any other sub compartment of a sensor cell. 
Such a prediction can clarify part of the noise generation. At the present time the noise in a 
conflned micro-domain cannot be directly measured. Instead, excised patch measurements 
are done fOl ^] , in which the cell structure is destroyed. Thus computation and simulation 
of mathematical models are the only tools for studying noise in this biological context. 



2 A Stochastic model of a non-Arrhenius chemical re- 
action 

2.1 Chemical reaction 

We consider two species of reactants, the mobile reactant M that diffuses in a bounded 
domain f2, and the stationary substrate S (e.g., a protein). The boundary dVt of the domain 
n is partitioned into an absorbing part dVLa (e.g., pumps, exchangers, another substrate that 
forms permanent bonds with M, and so on) and a reflecting part dVL^ (e.g., a cell membrane). 
In this model the volume of M is neglected. We assume that there are binding sites on the 
substrate. In terms of traditional chemical kinetics the binding of M to 5* follows the law 

M + Sfree^MS, (1) 

h 

where kf is the forward binding rate constant, kb is the backward binding rate constant, and 
Sfree is the uubouud substrate. 

However, when only a small number of reactant and substrate molecules are involved 
in the reaction, as is the case in a micro-domain in a biological cell, this reaction has to 
be described by a molecular model, rather than by concentrations. The description of this 
reaction on the molecular level begins with the following deflnitions: 

• M{t) = number of unbound M particles at time t 

• SA{x,t) = number of free sites in volume Ax at time t 

• S{t) = number of unbound binding sites at time t 

• MS{t) =M(0) - M{t) = number of bound M particles at time t. 

• So{x) = initial density of substrate 

• 'S'max = / Sq{x) dx = total number of binding sites in f2. 

Jq 
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The M(0) reactant particles are initially distributed with probability density mo (a;). The 
initial density So{x) integrates to S{0). We assume that the M particles diffuse and we 
denote by x{t) the random trajectory of an unbound M particle. We consider a small volume 
Ace about X, that contains initially so{x) Ax free binding sites and mo{x) Ax unbound M 
particles. 

The joint probability of an M trajectory and the number of bound sites in the volume 
Aa; is 

p{x, S,t\y)Ax = Pr {x{t) ex + Ax, SAix, t) = S \ x{0) = y} . (2) 

The function p{x, S,t\y) is the joint probability density to find an M particle and Sa{x, t) 
free binding sites at x at time t, conditioned by the initial position y of the M particle. 

The marginal probability density of an M trajectory is 

p{x,t\y) = p{x,S,t\y), 
free S 

where the sum is over all free binding sites in the volume Aa; . The number of free M 
molecules in the volume Aa; at x is 

Mf^gg(a;,t) = Mo Aa; ^ p{x, S,t\y)mo{y) dy. 

tree b 

The joint probability density function of x and S is 

p{x,S,t)= / p{x,S,t\y)mo{y) dy. 
Jn 

The two-dimensional process (Mfj,gg(a;, t), S{x, t)) is Markovian. The evolution of p{x, S, t), 
is governed by the diffusion of M particles in and out of the volume Aa; , and their binding 
and unbinding inside the volume Aa; . The influx in the time interval [t, t + At] is 
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J{x, S,t\y) ■ n{x) dSxAt = 

cJVolume Aa; 

/ [p{x,S,t + At\y) -p{x,S,t\y)]dx, 

^Volume Aa; 

which represents diffusion with coefficient D. Additional change in the contents of the volume 
Aa; is due to the binding and unbinding of M particles to the substrate. When there are 
S free binding sites in the volume Aa; (see Figure Q), the probability that one M particle 
binds to a free site in the volume Aa; in the time interval [t, t + At] is proportional both to 
5* and to the number M^^^^ix^t) of free M particles in the volume Aa; . It is given by 

kiMop{x,S,t)SAxAt. 

The probability that one particle unbinds in the volume Aa; in this time interval is propor- 
tional to the number so(cc)Aa; — S* of bound sites in the volume Aa; , given by 

A;_i[soAa; - S']At. 
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Thus the probabihty of S free sites when no change occurred in the number of free sites is 

p{x, S, t)Ax {1 - kiMop{x, S, t)SAx At - k^i [so{x)Ax - S]At} . 

The number of free sites can change to S at the end of the interval [t, t + At] if it was S + 1 
at the beginning and one bond was formed, or if it was S — 1 and one particle was unbound. 
The probability of this event, as calculated above, is 

fciMo Ax {S + l)pix, S+l,t)At + A;_i [so{x)Ax - S+ l)]At. 

The probability of a;, 5 at time t + At is therefore 

p{x,S,t + At)Ax = 

-V ■ J{x,S,t\y) Ax + 

p{x, S,t) Ax{l - MokiAxSp{x,S,t)At - k_i[soix)Ax - S]At} + 

Mok^iS + l)p\x, S + l,t) {AxYAt + 

k-i[so{x)Ax - S + l]p{x, S -1, t)AxAt. 
for S* = 0, 1, 2, ... , Smax, the coupled partial differential equations 



^M^^^A = -V-J{x,S,t)-K,p\x,S,t)S-k^,\So{x)-S]p{x,S,t) 
+ Ki{S + l)p'{x, S+l,t) + fc_i[S'o(a;) -S + l]p{x, 5 - 1, t). 



(3) 



where by definition J{x,S,t) is the joint probability flux at position x at time t, and S 
proteins are free. It is defined in the diffusion case by 

J{x,S,t) = -DVp{x,S,t). (4) 

The new (forward) binding rate is 

Ki = MokiAx 

and 

So{x) = sq{x) Ax = total number of binding sites in the volume Ax . 
Thus Ki is the probability flux into the binding sites. The boundary conditions on 5* are 

p{x, S,t) = for 5 < and S > So{x). (5) 

Remark 1. By summing equation Q over < S < So{x) and using the boundary condi- 
tions (0), we obtain 



for the marginal density of the M particles. It means that all bound and free M particles 
effectively diffuse. 

Remark 2. When at specific locations there can be at most one binding site, the system 
reduces to the coupled equations 

DAp{x, 0, t) - k-ip{x, 0, t) + Kip^{x, 1, t) 

(7) 

DAp{x, 1, t) + k^ip{x, 0, t) - Kip^{x, 1, t). 

Here Sq{x) can take the values or 1. When no molecules can escape from a bounded 
domain, the flux associated with p{x, S, t) satisfies the reflective boundary condition 

J-n\,^ = Q. (8) 



dp{x, 0, t) 
dt 

dt 



The initial condition, when no substrate is bound, is given by 0, 0) = mo(a;), hence 
p{x, 1, 0) = 0. When the total number of M particles stays constant (i.e., no particles leave 
the domain), adding equations ((Tj) gives in the steady state 

p(a;,0)+p(x,l)= (9) 



If the M particles can escape the compartment, e.g., by being absorbed in a part of the 
boundary dVLa, the condition (jSJ should be changed to 

J ■ ^lan-an^ = (10) 

and 

pix,S,t)\,^^ = 0. (11) 

In this case (0) no longer holds. 



Remark 3: Obviously, So{x) takes only integer values. We assume that its discontinuities 
are located on smooth interfaces. The density p{x, S, t) and the normal component of the 
flux J{x, S, t) ■ n{x) are continuous across the interfaces for all S* = 0, 1, . . . , S'max- 



2.2 Moments of the pdf 

Statistical moments of the pdf can be computed from equation The average and the 
standard deviation of the number of bound proteins are evaluated for equation ((Zj). The 
mean of the number of bound proteins at time t is given by 

{Sh{t)) = [ Soix)pix,l,t)dx. (12) 
Jq, 
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The standard deviation is given by 

a\t) = (Slit)) - {S,{t)y = j^Sl{x)p{xXt)dx- (^J^So{x)p{x,l,t)dx^ . (13) 

When the proteins are uniformly distributed on a subset Q' C Q, the distribution So{x) is 
given by the characteristic function of the subset Q', 

So{x)=xn'{x)j^^, (14) 

where the total number of binding sites is 

No= So{x)dx. 
Jn 

We obtain the standard deviation of bound sites from the expression fll3|l as 

^'(t) = {sm-{s»{t))' (15) 

^0 ^2,.,\ f No^^ 



Soit)-So\t)) = (-^] 5o(t)(l-So(t)) 



where 



S'o(t) = / p{x,l,t)dx (16) 
Jn' 

is the fraction of bound sites. Note that 

So{t) < 1. 

2.3 Standard deviation of the number of bound protein in some 
cases 

We consider the one-dimensional case where fl = [0, L] and 5*0(0;) is either or 1 in intervals. 
In the steady state the system ((7j) is 

= DAp{x,0)-k_iSo{x)p{x,0) + Kip^{x,l) 
= DAp{x,l) + k^iSo{x)p{x,0)~Kip^{x,l) (17) 
^ = P{x,0) +p{x,l), 



which reduces to 



Dp"{x,l) + k_^iSo{x) 
We convert to densities by setting 



--p{x, 1) 



cm{x, 1) = Mop(x, 1) 



Then fllHj) becomes 

Dc'j^{x, 1) + A;_iS'o(x) 



- CAf(x, 1) 



Mn 



c^^(x, 1) = 0. 



(19) 



The function cm{x, 1) is supported where the protein are located. 

When the boundary conditions are reflective for the M trajectories, using the uniqueness 
of the solution, then at a point x, where 5*0(0;) is supported is 



cm{x, 1) 



2A;_i5'o(a;)Mo/L 



k-iSo{x) + v/(A:-i5o(x))2 + 4iriA;_iS'o(x)/L' 



In particular. If the substrate is uniformly distributed in intervals Sq{x) = Ns/L and Mq ^ 
Ns, then the fraction of bound M particles is 



^^0 Jo 



1 + a/1 + 4 



MokiAx 



1, 



(20) 



which means that practically all M particles are bound. If Mq ^ Ns, then gives 

-L 



Jo 



Moki Ax 



< 1. 



(21) 



In this case, the variance of the fluctuations in the number of bound M particles, which 
is the same as the number of bound sites, as a function of Mq and Ns, is given by 



cr|(Mo) =Pm(1 -Pm) 



1 + 4 



MokiAx 
Nsk.i V 



1 - 



1 + 4 



MokiAx 
Nsk^i / 



(22) 



The graph of (t|(Mo) vanishes for both Mq and Mo 00 and has a unique finite 
maximum, as in figured 



8 



2.4 General equations when binding proteins are located on the 
boundary of the domain 

We consider the same problem, but with binding sites located on the boundary dQ with 
surface density 5*0(3;) (see figure El). This may represent for example binding to gated chan- 
nels, which open when they bind agonist molecules. Some channels may need to bind several 
agonist molecules to open, however, we consider here the case that a single agonist molecules 
opens the channel upon binding. We assume that the agonist molecule cannot escape fl. 

The initial Mq agonist molecules diffuse in Q and are reflected at dQ at non binding 
sites, but can bind to a free binding site on a protein channel in the membrane with a 
certain forward binding rate ki. When this occurs the channel opens and stays open as long 
as the agonist is bound. The bound agonist is released from the bound state at a backward 
rate k_i. Then the number of open channels is the number of missing M particles in Q. We 
write 

Mo = / c{x,t)dx+ / S^Q^^^{x,t)dSx. 
Jn Jan 

The equations for the density of M particles in Q and in dQ is 

^ "^bound^^' ^) 

XedQ 

= kic{x,t)[So{x) - 5bound(*'^)] ' ^-i'^bound(*' (23) 
The probability of k bound sites on the boundary at time t satisfles the equations 

Pk{t) = -Pk{t)Ki(f c{x,t)[So{x)-S^^^^^{x,t)] dSx-kPk{t)k.i 
Jan 

+Pk-i{t)K^ (f cix,t) [So{x) - 5bound(^'^)] ^Sx + (k + l)Pfc+i(t)A;_i, ,(24) 
Jan 

Po(t) = -Po{t)K,(fc{x,t)[So{x)-S^^^^^{x,t)]dSx + Pi{t)k_u, (25) 
Jan 

Ps,it) = -Psoit)K^ (f cix,t)[So{x)-S-^^^^^ix,t)] dSx-SoPso{t)k_^ 

Jan 

+Ps,-l{t)K^ (f c{x,t) [So{x) - 5bound(^'^)] (26) 
Jan 

where, according to eq. (f^ 

'^bound(*'^) = kic{x,t)[So{x) - 5bound(^'^)] ~ ^-i'^bound(^' 



-D 



dc{x, t) 
dn 



J{x,t) 



nix] 
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Here < k < So, where 

So= f So{x) dSx, 
Jan 

that is, P_i(t) = Ps„+i(t) = 0. 

The moments of the number of bound sites are 

So So 
k=l k=l 

The variance in the number of bound sites is 

At) = {Slit)) - {s,{t)r. 

3 The fluctuation in the number of particles in a push- 
pull chemical reaction 

A push-pull chemical reaction consists of a source that produces particles at a given rate and 
a sink or a killing term that destroy or remove the particles from the system with its rate. 
When the movement of the particles is driven by diffusion and the source and the sink are 
not uniformly distributed, the number of particles fluctuates. In this section, we propose an 
approach to estimate the mean number and the fluctuations. A permanent regime imposes 
that the rate of the sink and the source satisfy some specific conditions. At equilibrium, in 
the limit of a large number particles, the push-pull mechanism reaches a steady state and the 
number of particles does not fluctuate and can be computed using the rate constant. But 
for a small number, the analysis requires to study separately the dynamics of the particles 
and especially the law of injection by the source. 

In the neurobiological context, many biochemical reactions are based on push-pull mech- 
anisms: at a given location, a source produces molecules and somewhere else, an enzyme 
modeled as a sink, destroys the molecules with a certain efficassy. In general, sinks and 
sources are not uniformly distributed, which induces fluctuations. As an example, gating 
molecules that can open channels are produced by such push-pull mechanism, and the fluc- 
tuations in the number of molecules induce a fluctuation in the number of open channels, 
which is reveals at the cellular level. It is of particular interest to examine the situation 
where molecules move by diffusion, the enzyme sink is represented by a single molecule and 
the sources are uniformly distributed. For example, this occurs inside a compartment of 
photo-receptors for the regulation of cGMP molecules. 

3.1 The push-pull mechanism 

The sink of the push-pull mechanism is modelled by a killing term located at a single point. 
The killing term is a Dirac function, ( see ^T] for the exact mathematical interpretation). 
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The source is assumed to be uniformly distributed and molecules are produced at a constant 
rate 7. 

We assume that the molecules are independent and their movement in a domain Q can 
be described by the stochastic differential equations 

Xk = b{xk) + V2Dwk, k = l,...,N, (27) 

where b{x) is a drift vector. Since the molecules cannot escape, reflecting boundary condi- 
tions are imposed at the boundary dfl. If a particle is injected in Q at time t = 0, its pdf 
p{x, t) satisfies the Fokker-Planck equation 

dp{x, t) 



dt 



J ■ n\ 



-V • J{x, t) — ki5{x — Xi)p{x, t) 



{2i 



an 



0, 



where the probability flux density vector is given by 

J{x, t) = -DVp{x, t) + b{x)p{x, t). (29) 

At any moment of time, when there are particles inside Q, the concentration c{x, t) is 
given by 

c{x,t) = Np{x,t). (30) 

Particles are injected with a Poisson stream at a rate 7, that is, at independent and identi- 
cally distributed inter-injection times, whose pdf is 

The mean inter-injection time is E{T) = 7^^. The probability that a molecule injected at 
time t = at a point y survives at time t is 

Sy{t) = Pr{x{t) eQ} = p{x,t\y)dx. (31) 

Jn 

When y = 0, we denote S{t) = SQ{t). To compute the mean and the variance of the number 
of molecules N{t) surviving in Q at time t, we use here the renewal equation |H 131 1^. 



PT{N{t) = 0} 
Pr{A^(t) = 1} 



PT{x{t) ^ Q} ■ 
FT{x{t) G n} ■ 



f{s) Pi{N{t-s) = 0}ds + 



FT{N{t) = n} 



f{s)FT{N{t-s) = 0}ds + 

+ PT{x{t) iVt}- f f{s) Vi{N{t - s) = 1} ds. 
Jo 

PT{x{t) en} - [ f{s) PT{N{t -s)=n-l}ds 
Jo 

+ PT{x{t) ^Q} ■ [ f{s)PT{N{t-s) =n}ds, n > 1. 
Jo 



f{s)ds 
f{s)ds 



(32) 



(33) 



(34) 
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The expected number EN{t) of molecules surviving in Q at time t is 

oo 

EN{t) = ^nPT{N{t) = n} 

n=l 

(ft POO 
J fLis)EN{t- s)ds + J f{s)ds 

+ Pi{x{t) ^Q} ■ [ f{s)EN{t- s)ds 
Jo 

= PT{x{t)eQ}+ [ f{s)EN{t- s)ds. 
Jo 

The integral equation (j35|l is solved by the Laplace transform as 



(35) 



n(r) ^ (36) 

where S{t) is the Laplace transform of Pr{a;(t) G Q} when the initial position of insertion 
is 0. 

fir) = (37) 
7 + r 

and S{t) is the Laplace transform of the survival probability. Therefore the Laplace trans- 
form of EN{t) is given by 

nir) = tmM. (38) 



3.1.1 Computation of S{r) in a driftless one-dimensional model 

We consider the one-dimensional equation ()28|) in f2 = [0, L] and b{x) = 0. We assume, for 
simplicity, that L = n and D = 1. 

PL 

The survival probability given by S{t) = / p{x,t \ y)dx is also (see |llj) equal to 

Jo 

S{t) = l-ki [ p{x,t\y)dx. (39) 
Jo 

The Laplace transform is 

Sir) = ^-f^^Pi^^'-\y)_ (40) 

T 

p{xi, T I y) can be computed using the Green's function for the Neumann problem for 



given by 



2 ~ 



G(x,t\y) = l-\ — cosnx cosny. (41) 
7r ^-^ 



n=l 
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Following the computations of JT] , an integral representation of the solution of equation | 
is given by 



p{x, t\y) = —ki / p{xi, s I y)G{x, t — s \ xi) ds + G{x, t\y). 
Jo 

The Laplace transform of equation ()42|) is given by 

p{x, T\y) = -kip{xi, T I y)G{x, r \ xi) ds + G{x, T\y). 

Thus, 

G{xi,T\y) 



p{xi,r\y) 



1 + kiG{xi,T I xi) 



and 



With 



' ki 1 + kiG{xi, T I Xi) 



1 2 cos nx cos n?/ 



G{x,r\y) = - + - V 
r TT — ^ 



n=l 



?T,2 + r 



for x,y G (0, vr), and to first order in r, we obtain 



-E 

1 ^ _ 
p(x,r|?/) = — + 
fci 



cos no; cos n?/ cos nx 



, 2ki cos rix cos ny 
ki + > - 



+ o(r) 



n=l 



1 



27r^ ^ ^ ^ 



, , 2 /vr^ TT x2 

l + - xH 

' TT V 6 2 2 



r + o(r). 



Hence 



1 - kip{xi,T\y) 



(42) 



(43) 



(44) 



(45) 



(46) 



271 



{x' - y') + fcr^ 



2 /tT^ TT X2 

a; H 

TT V 6 2 2 



+ 0(1). 



Using the normalization condition that 5(0) = 1, we find the long time asymptotics 

5'(t) ~r^+oo exp(-Q;t), (47) 
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where 



, 2 /tt^ tt 

1 + - x-\ 

TT V 6 2 2 



The un-normahzed time constant for x,y & [0, L] is 



-TT 
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3.1.2 The steady state limit 

In the steady state, the mean number of molecules N{t) surviving in at time t is 



hm EN{t) = lim Tri(T) = 7^(0) 



Note that 



poo nL 

S{0) = / p{x,t\0)dx. 
Jo Jo 

If at time t — the particle is injected randomly at a point y, we index S hy y and 

Sy{^)- / / p{x,t\y)dx^ / G{x\y)dy 
Jo Jo Jo 

where G{x \ y) is the solution of 

-5{x-y) = D?^-k5{x-x^)G 

dG{Q I y) ^ dG{L I y) _ ^ 
dx dx 

Thus, 

where 9{x) is the integral of the Heaviside unit step function. It follows that 
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At steady state, the mean number of molecules surviving in Q at time t is given by 

iV,M^,5,(0)^.(-(^g^£+(^ + ^) (51) 

dx 

If the particle is initially uniformly distributed, Pr{y e [x, x + dx]} — — , then the steady 
state mean number of surviving molecules is given by 

iV(co)=,£,[S,(0)l=,(-^+(^ + Q. (52) 
3.1.3 Variance 

The second moment of the total number of surviving particles is defined as 

oo 

EN^{t) = J^n2pr{Ar(i) = n}. 

n=l 

We will now compute such number using the renewal equations. We have 

EN\t) = Fr{x{t) en} [ f{s)k'^ Fr{N{t - s) = k - 1} ds 

Jo 

+ Pr{x{t) i [0, L]} / f{s)k^ VT{N{t - s) ^ k} ds + Fr{x{t) ^Vt} j f{s) ds, 
Jo Jt 



which leads to 



J poo pt 
' f{s)ds + Fi{x{t)en} j EN^{t- s)f{s)ds + 
t Jo 

PT{x{t)^[0,L]} [ EN'^{t- s)f{s)ds 
Jo 

+ J2 P^{x{t) e ^} / (2A; + l)Pr{N{t) = A; - ds. 

Thus, 

/t POO 
EN^{t - s)f{s) ds + Fr{x{t) e il} f{s) ds + 

2 Pr{x{t) eQ} [ E{N{t - s)f{s) ds + Pr{x{t) e Q} [ f{s) ds 
Jo Jo 



and 

ft 



EN'^{t) = [ EN^{t - s)f{s) ds + PT{x{t) e Q} + 2 Pr{x(i) e n} [ f{s)E{N{t - s))ds.{53) 
Jo Jo 



15 



Remark. 

The variance can also be written in the form 



= EN\t) - \E{N{t))Y = Fr{x{t) e Vl} - Fi{x{t) e ny + 



EN^{t - s)f{s) ds - ( [ E{N{t - s))f{s) ds 
\Jo 



3.1.4 The steady state variance 

To compute the variance in the steady state hmit, as t — > cxd, we Laplace transform equation 
(ESI), 



It follows that 



S{r) 



1 - /(r) 1 - /(r) 7o 



S{t) I f{t-s)n{s)ds 





e"'" dt. 



EN^ (oo) = lim EN^{t) = lim™2(r). 

t—*oo r— >0 

Writing 

Tn^{T)=I{T)+II{T), 

we obtain the short r approximation 

S{r) 



It) = T 



1 - m 



To evaluate 



II(t) = T- 



1 - fir) Jo 



S{t) / f{t-s)n{s)ds 



e"'^ dt, 



we recall that f{t) = "'^ and approximating S{t) by e and interchanging the order of 
integration, we obtain 



//(r) = r7 



c\ nco POO 

1 - fi-r) Jo Js 



r7 1— = / n(s)e-*("+")cit, 

^l-/(r)7 + r + «yo 

ri- TTT — : : — "^v + ^) 



1 - /(r) 7 + r + a 
7 + a 
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Finally, 

1 



and 



EN'(oo) = 7^(0) + 2Y n(a) 

7 + a 



= EN\oo) - {EN{oo)f = 7^(0) - 72^^(0) + -^^n{a), (54) 

7 + a 

where S{0) = — i^^^ the time a depends only on ki, the diffusion constant and the length 

7 

7 -|- G; - - 

L, as given in formula (jlHI), and n(a) = S{a), where S{a) is the Laplace transform of 

a 

the survival probability at time a. 



3.2 Push-pull chemical reaction in the continuum limit 

When the number of reacting molecules is large enough, regulated by a push-pull mecha- 
nism of hydrolysis and synthesis, the reaction-diffusion equations is sufficient to describe the 
evolution of such system. We assume that the molecules diffuse inside a domain Q, but the 
synthesis occurs only uniformly in a domain Q' G Q while the hydrolysis is performed by 
isolated enzymes. The hydrolysis can be modelled as a killing measure k{x) (see When 
there are a discreet number Nh of killing sources, the killing measure k{x, t) is the sum of 6 
measures given by 

k{x) = "^kidix - Xk), 

1 

where ki is the arrival rate constant to the killing sources. If the boundary of the domain is 
divided into an absorbing part dQa and a reflective part dflb, the concentration c{x,t) can 
be defined as the probability density function p{x, t) of the system of N particles as 

c{x, t) = Np{x, t). 



DAc{x, t) - k{x, t)c{x, t) + 7Xq' 

0, 

0, 

Co(x), 



and 



dc{x, t) 
dt 



dc{x, t) 
dn 

c{x, t) 



c{x,0) 
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where 



/ CQ{x)dx = N. 
Jn 

The mean number of particles N{t) at any moment of time is given by 



N{t) = / c{x,t)dx. 

When the particles are contained inside a fixed domain and cannot escape, is empty. The 
steady state density c{x) = lim^^oo c(a;, t) is a solution of 







DAc{x) - k{x)c{x) + 7Xn' 
dc{x) 



dn 



In one dimension Q' = Q = [0,L] and the killing measure is reduced to a single point. The 
steady state equation is 











^~dx^ ^^^^^ ~ a;i)c(x) + 7, 

dc{x) dc{x) 



(55) 



dx 



x=Q 



dx 



x=L 



where ki is the forward rate for one active site ( it has units of length per time), 7 is the 
uniform injection rate of particle due to synthesis ( it is the number of injected particles per 
second per unit length). Two integrations of equation lead to the solution 



c{x) = 7 



L6{x — Xi) 
D 



X 

2D 



(-- 






2Dy 



(56) 



where 6 is the integral of the Heaviside function: 

{0 for X < 
X for X > 

Finally, the total number of free molecules is given by 



which should be compare to the steady state mean number obtained in the discrete case 
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4 Markovian model of cell noise 



We propose in this section an alternative model for computing the fluctuation of the number 
of bound molecules, based on a Markovian approach. We consider a domain f2 C with 
M mobile agonist molecules and S receptors, embedded in the boundary dVL, that open 
channels when they bind a single mobile agonist molecule. We assume that the receptors 
occupy a small portion of the surface area of dVL. The agonists diffuse in VL independently of 
each other. Bound agonists are released independently of each other at exponential waiting 
times with rate 

For a single receptor and a single agonist the time to binding is the first passage time 
to diffuse to a small portion of the boundary, dVta-, which is absorbing and represents the 
active surface of the receptor, whereas the remaining part of dVt is reflecting. It can be 
shown that the probability distribution of the first passage time to dVta is approximately 
exponential with rate 



1 



in) 

where (ri) is the mean first passage time to dQa- 

When there are S channels, k{t) of which are free at time t, and M agonists (gating 
molecules) , 

N{t) = {M - S + k{t))+ (57) 
of which are free to diffuse in Q at time t, where 

= max{0, x}. 

Thus 

{S - M)+ < k{t) < S. 

We assume that the pdf of the time for the next receptor to bind is well approximated 
by the exponential pdf with instantaneous rate 

Nmt) 

This assumption is justified if the total area of the absorbing boundary (the channels) is 
small relative to the surface area of the reflecting boundary PJ. 

It follows that when N{t) = N and k{t) = k, the mean time to bind is 

Nk k{M-S + k)+ 



in) in) 



The number of bound receptors at time t is a birth-death process with states 0,1,2,..., min{M, 
and transition rates 
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The boundary conditions are 

^s-*s+i = 0, Ao^-i = 0. 

Setting 

Pkit) = PT{k{t) = k}, 

we have [22] 

P(5_A/)+(t) = -k^lSP(^S-M)+it) + \iP(S-M)++l{t) 

Pk{t) = -[\k + k^i{S~k)]Pkit) + \k+iPk+iit) + k^i{S-k + l)Pk^,{t) 
for k = {S- M)+ + 1, . . . , ^ - 1 

Ps{t) = -XsPs{t) + k^iPs-i{t). 

The initial condition is Pk,q{0) = Sk^s^q,o- In the steady state the average number of open 
channels is 

5 

(koo) = jPjy 

j=(5-A/)+ 

where Pj = lim Pj(t), and the stationary variance in the number of open channels is deter- 

t — >oo 

mined from the second moment 

j=(S-M)+ 

by 

a'{M,S) = {kl)-{k^f. 

4.1 The steady state approximation 

In the steady state 

= —k^iSP(^s-M)+ + ^lP{S-M)++l 

= -[\k + k^i{S-k)]Pk + \k+iPk+i + k^iiS-k + l)Pk-i 

for k = {S -My + 1,...,S -1 
= -XsPs + k^,Ps-i, 

20 



which gives for {S — M)~^ ^ k < S 



p nt5-.+i A. _ p nt^-.+i ^iM-s+ 1)- 

S-k — 777% — ^S- 



„ _ „ Y[^i={S-M) ++l^i _ „ nf=(S-M)++l ~ "S* + 



(5 _ (5 _ M)+)\k'_'^^'-''^" '{S-{S- My)\{k^M)Y 

The constant Pq is determined from the normahzation condition 

s 

k={S-M)+ k=S-(S-M)+ 

Thus, 

- ^ 



^ ^k=l u\Uk i + / TTT^ — 7-; 



^'^-1 4^ fc!((ri)A;_i)'= 



k=l 

{S-M)+ {S-M)+ y-fS 



(kl) 



k=S-l k=S-l ^' ^' 

k=S-l k=S-l '^^ ^' 



4iM) = {kl)-{k^)\ (5^ 
The graph of (7|(M) vs M is given in Figure El 



4.2 The steady state mean number of open channels 

Using the above Markov model an exphcit formula can be obtained for the mean number of 
open channels ( that bind an agonist molecule). It is assumed here that a channel can bind 
only a single agonist molecule at a time. This approach does not give any information about 
the fluctuations. 

The steady state fraction n\ of bound agonist molecules in a bounded domain Vt that 
contains S channels and M gating agonists is 

(60) 
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where ki is the forward binding rate in solution and k_i is the backward rate. If Sch is the 
total surface occupied by the channels, then Sch = SS^f^, where S"^^ is the effective surface 
occupied by a single channel. 5*^;^ is defined by the area where the electrostatic force of the 
channel is sufficient to bind an agonist. Recall that ki is by definition the mean number of 
agonist molecules arriving at the channel per unit of time. Using the result of the previous 
section, we can write 

Thus, using equations ()60p and ()62j). we find that the mean number of bound agonist, which 
is the same as that of bound channels, is given by 

nl = -J M, (63) 

l + ^{M-nl){S-nl) 

Och 

with 

where \Q\ is the volume of the Q and \dQ\ is the surface area of its boundary (see Pj). For 
nl small relative to M and S, 

^ M (65) 



Tik 



-MS 



Sch 

(66) 

^ M, (67) 



"^ch 



which gives explicitly the number of bound channels as a function of the geometrical pa- 
rameters. A similar formula can be derived, when a single channel can be bound to several 
gating molecules. 



4.3 The Michaelis-Menten law in micro-structures 

The rate of production of a product P from a substrate M by a catalytic enzyme E is usually 
described in text books by the Michaelis-Menten law for reactions in solution. In confined 
micro-domains, the above analysis can be used to estimate the number of P produced in the 
reaction. 

The chemical reaction is described by 

kfi 

M + E ^ ME ^-^ E + P. (68) 

hi 
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and a master equation for the joint probability that the number of P molecules produced is 
k and q enzymes are bound at time t, Pk,q(t) = Pr{P(t) = k, E(t) = q}, can be derived as 
above. In particular, the kinetics of the reaction (|68p are 

Pk,q{t) = (1 - q{k-l + k2) - XM,q)Pk,q{t) + H^l + I) Pk-l,q+l{t) + k.l{q + l)Pk,q+l{t) 
+ XM,q-lPk,q-l(t), ioT 1 < q,k < Eq - 1 

P0,q{t) = (1 - XM,q)P0,q{t) + ^-M + l)^0,,+lW + AA/,q-lPo,,-l (t) , 
Pk,Eo{t) = (1 - k{k-l + k2))Pk,Eo{t) + k2qPk-l,Eoit) + XM,Eo-lPk,q-lit), 

Pmo.oW = (1 - Mo(fc_i + k2))PMo,oii) + hPMo-ifiii)^ 
with the initial condition Pk,q{0) = 5m-o,o and 

_ M{Eo-k) 

^M,q — 5 

n 

M = {Mo-k-q) + . 

Here Mq is the initial number of substrate molecules, M is the number of unbound substrate 
molecules, and Eq is the total number of enzyme available. The mean number of P produced 
at time t is 

{P{t)) = J2 kPr{P{t) = k} = J2kJ2 ^^{^(^) = ^' ^ W = (69) 

k k q 

The only steady state solution of equation (j69j) is zero as any steady state. 

4.4 Fluctuations in a push-pull system with binding 

We consider now a push-pull system, where gating molecules can also bind and unbind to 
some proteins. We consider two main approaches to the description of this dynamics. The 
first approach is based on the procedure used in the first model in Section |21 and consists in 
deriving an equation for the joint probability of a trajectory, the number of bound sites in 
the volume Ax, and the total number of molecules in Q, 

p{x, S,M,t\ y) Ax = Pr {x{t) E x + Ax, Sa{x, t) = S , M\ x{0) = y} . (70) 

The function p{x, S, M, t\y) is the joint probability density to find an agonist and S free 
binding sites at cc + Ax at time t, and M agonists in the domain, conditioned by the initial 
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position y of the agonist. An analysis similar to that of Section |21 leads to 
dp{x,S^M,t) ^ . j[x,S,M,t) - K^Mp\x,S,M,t)S - k_i[So{x) - S]p{x,S,M,t) 

+ KiM{S + l)p\x, S + l,M,t) + A;_i [So{x) -S + l]p{x, S-l,M,t) + 

-fp{x, S, M, t)-{M - SyK^ip{x, S, M, t), (71) 

where by definition J{x, S, M, t) is the joint probability flux at position x at time t, and S 
proteins are free for M agonists. It is defined in the diffusion case by 

J{x, S, M, t) = -DVp{x, S, M, t). (72) 

The new (forward) binding rate is 

Ki = kiAx, 

7 is the uniform production rate, and K_i is the uniform killing rate for the (M — 5')"'" free 
agonist molecules. The moments can be calculated from p{x, S, M, t), as above. 

The second procedure consists in adding directly the push-pull effect in the Markovian 
model. Defining the joint probability that k channels are free and M{t) are in the reacting 
domain at time t, 

Pk,,{t) = PT{k{t) = k,M{t)=q}, (73) 

neglecting the distribution of the sources and sinks responsible for the fluctuation of the 
number of molecules, we consider only the case that hydrolysis and synthesis occur at 
exponential waiting times (i.e., are Poissonian), with rates K_i and 7, respectively. Following 
the same steps as in Section HI the master equation for Pk,q{t) becomes 

P{S-M)+,q{t) = -k-lSP^S-M)+,q{t) + \l,q{q - S+ 1)^ P^S- M) + + l,q{t) 
Pk,q{t) = -{Xk,q+{S -k)k_l+-i + K_iq)Pk,q{t) + 

Xk + l,gPk+l{t)) + k^l{S -k+ l)Pk-l,g{t) + -fPk,q~l + K.l{q + l)Pfc,5 + l, 

for A; = (^-M)+ + l,...,^-l, 
where we recall that 

Nk k{M-q + k)+ 



X 



k,q 



in) in) 

The initial condition is Pk,q{0) = Sk^s^qfl. 
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Finally, if the push-pull rate is much slower than the binding rate, the previous results of 
the Markovian model can be used directly and the mean number of free channels, N^{t), 
is directly computed using the moments and 

NUt) = 5^(fcoo)(M) PT{N{t) = M}, 

M 

where {koo){M) is the mean number of bound channels when there are M gating molecules 
in the domain, 

(S-M)+ 

(k^){M)= iS-kyPs-k{q)- 

k=S-l 

Here, by definition, the dependence of Ps-k on q is that obtained in equation (jSHI), and 
Pr{A^(t) = m} is the probability that M molecules are in the reacting domain at time t (see 
Section In the limit t ^ oo 

(iV^M) = J](A;„o)(M)Pr{iV(oo) = M} 

M 

(iV^(oo)) = 5^(A;o,)2(M)Pr{iV(oo)=M}. 

M 



5 Conclusion and biological implications 

5.1 Comparison of the models 

We presented here two models that describes fluctuation due to binding and binding of 
agonist to some fixed proteins. In the Markovian model, the only geometric feature of the 
cell that enters the model is the cell's volume. The distribution of the channels, as well as 
other geometric features are ignored. The advantage of such approach is that it gives explicit 
estimates of the mean and the variance (see equation I58|). In the first model, more details of 
the geometry and organization of the channel are captured, however, the computation of the 
moments requires the solution of a system of partial differential equations. The expression 
of the variance is given in general by equation [121 and in steady state by expression |221 



5.2 The Forward binding rate 

The Markovian model does not rely on the forward binding rate and it is replaced here by 
the mean arrival time to the channels (ri). The initial forward binding rate per molecule is 
the arrival rate of a gating molecule to any one of the binding sites on the proteins, and 
can be defined as 
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The number can be, for example, the number of proteins in a given volume. It can be 
computed from the concentration [A^]. In that case Ki can be rewritten as 



The Markovian approach proves that the traditional forward binding rate has to be aban- 
doned and a dynamical rate has to be used instead. The first model uses another definition 
of the forward binding rate, that can be related to the traditional one. 



5.3 Biological implications 

The mean and the variance of the number of bound molecules were derived in the first model 
for a finite system of molecules. In particular formula (fT3j) proves that the variance depends 
on the total number of bound molecules at time t So{t) ( see formula (fTB|) ). Note that 5*0 (t) 
is the integral of the joint probability density function p{x, l,t), which is a solution of a 
partial differential equation. Thus the number of bound molecules is a complicated function 
of the geometry, the distribution of the substrate molecule and the interaction between the 
molecules through the number of binding sites. 

As seen in the experimental data [121 HH CZl UH] , the variance (fT3)) does not only reflect 
binding to proteins, but also a dynamical process of binding and unbinding in a microstruc- 
ture, which involves the geometry, the distribution of the protein and the diffusion process. 
Only at high concentration of binding molecules, a molecule that leaves a binding site will 
be immediately replaced by another gating molecule. In the regime of high concentration, 
when the protein represents a gating channel located on the cell membrane, the fluctuation 
of the current represents an intrinsic property of the channel. The fluctuations, then, are 
proportional to the gating property of the channel. But for many sensor cells, the concen- 
tration of gating molecules is not high, the small number of gating molecules is the cause of 
perpetual fluctuations due to binding and unbinding. 

The graphs 1^ and El of the current noise variance vs So{t) has almost the shape of an 
inverted parabola, which implies that there it has unique maximum point at a specific 
location. In the steady state limit, it is obtained in case of the first model for SQ{t) = 1/2. 
Biological microsystems, such as photo-receptors, usually operate at low noise levels, far 
from the maximum of the graph. 

Controlling the membrane voltage fluctuation for sensor cells is a crucial issue for the 
transduction process that consist in transmitting a molecular signal. Channels fluctuation is 
not responsible only for the noise of the cell, but the biochemistry underneath and the orga- 
nization of the cell play a crucial role. In a set of experiments P^ITHllT7| ll9j. the membrane 
noise was measured but the purpose was to identify the properties of the channel rather 
than the molecular dynamics. Using the results of the previous mathematical analysis, it is 
a hard problem to convert the data of fluctuation obtained for a detached patch experiment 
to the fluctuation inside a single cell. A separate mathematical analysis has to be performed 
to identify such fluctuation. 
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In the experimental data, the variance of the fluctuation has been related JH] to the total 
current, when the probability of the channel to be in an open state is independent from the 
other. When the gating molecule is not as abundant, channels share such resource and when 
a gating molecule is bound to a channel, this resource is not available for the neighboring 
channels. This effect coupled the probability of the channel to be in open state through 
the number of gating molecules. It is useful to compare the experimental result El with the 
simulations of the theoretical model, obtained in figures ^ and |21 This comparison suggests 
that the tail distribution of the variance cannot be approximated everywhere by a parabola, 
as it was done, under the assumption that the opening of the channels is independent. This 
assumption breaks at high concentration, where the channels are coupled through the gating 
molecules. 

In a system composed by a cell membrane containing channels and gating molecules only, 
two time scales should be dominant. The first time scale is related to the time it takes for 
a gating molecule to find a channel. This time depends on the geometry of the domain, the 
number of gating molecules and the number of channels: for few molecules, this time can 
be approximated by the mean time it takes for a gating molecule to find a small absorbing 
boundary. This problem has been treated in jH] and this approximation was used in the 
Markovian model of the paper to estimate the number of open channels. This assumption is 
valid because channels or binding proteins occupy a small portion of the boundary surface. 
The second time scale is related to the backward binding constant of the gating molecule to 
the channel. The backward binding constant depends on the property of the channel only 
and thus does not depend on the statistical property of the system. 
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Figure 1: Chemical reaction inside a microdomain. Each free particle (white circle) moves 

according to a Brownian motion until it binds to a free site (black color). An occupied site cannot 
accept any other particles. A molecule is composed of a finite number of sites. Parts of the 
microdomain is reflective and the rest is absorbing. 
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Figure 2: Schematic reaction in a microdomains where channels are located on the 
boundary. Molecules move according to a Brownian motion inside the domain. When a particle 
hits a free channel, it can bound and stays there for a certain amount of time. Particles can escape 
through the absorbing boundary. 
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Figure 3: The variance square equal to o"|(M), for various values of the backward binding rate 
A;_i = [250,500,1000,1500]. The total number of substrate molecules is fixed at 5 = 10 and 
T = 0.01ms. 
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Figure 5: Experimental analysis of the fluctuation in the number of open chan- 
nel. The left panel corresponds to the current response for various concentration of the gating 
molecules(cGMP), in a detached patch experiment of a cone Photoreceptors membrane. On the 
right panel, the variance is plotted versus the mean current, which is proportional to the number 
of gating molecules (cGMP). This picture was published in ^Hl figure 2. 
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